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Abstract 

Exact analytical expressions in terms of generalized confluent hypergeometric func- 
tions for the transition amplitudes of neutrino oscillations in presence of matter are 
computed for an arbitrary number of species. The density of matter is assumed to be 
exponentially decaying. The results can be used for the description of matter-induced 
| neutrino oscillations in the Sun which can take place when the solar neutrinos prop- 

agate radially from the interior to the surface. Expressions are particularly simple in 
| the limit of infinite propagation time as is suitable for the case of detection at Earth. 

^ ' PACS: 14.60.Pq, 02.30.Gp, 02.30.Hq 
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1 Introduction 



The explanation of the solar neutrino flux problem by maximally mixed neutrino oscilla- 
tions in vacuum or by resonance-enhanced oscillations without large mixing (MSW effect) 
is well known by now and it has been treated abundantly in the literature ([|l], [|, |||). 

For stable and relativistic neutrinos v = (ui,i = 1,N U ) of equal momentum the prop- 
agation in vacuum is described in the basis of mass eigenstates by 

%d t v = H°u 

where 

H° = DiagiEi = J{p 2 + m 2 ) ~ p + m 2 /2E) (1) 

The common constant p can be absorbed as an unimportant phase and is usually omitted. 
The weak eigenstates v' are related to mass eigenstates by (it* for antineutrinos) 

V = U V 

The transitions between v' are described by the amplitude matrix 

A = ue~ m0t u~ l 



In presence of standard matter with arbitrary electron number density the propagation 
is usually well aproximated by 

idtv = {H° + P {t)uAu- l )v (2) 

where A is a matrix with An as only non-zero element. p(t), essentially a forward scat- 
tering amplitude, is proportional to the electron number density of the medium 

p(t) = ±y/2GpN e (t) 

The plus/minus sign is for neutrinos/antineutrinos respectively. 

The difference of the eigenvalues of H° ( the inverse of the usually defined oscillation 
length) is typically considered to be, in the solar case, of the order 

mf-m 2 0.1 -leF 2 fi 7 

- - 3 w — — — — ps 10~ 6 - 1(T 7 eV 

2E 1 MeV 

The electron number density can be parametrized (data taken from H ) , for sufficiently 
far distances from the solar core as 

N e (r) = iV exp(-Ar) = 160 ± 30exp(-10.6 ± 0.2 r/r ) mol/cm 3 

with r, ro the distance from the center and solar radius respectively. At the solar center 

N e (0) ps 100 mol/cm 3 
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So the value for p(t) varies between zero and its maximum value of ~ 10 -11 eV at the 
center. 

Later we will see that the following dimensionless quantities play an important role: 

1 3 1.75 10 8 — i- J (« 1 for E= 1 MeV and Am 2 = 1(T 7 eV 2 ) 
2EX E (Mev) v ; 

p(t)/X varies between 5.5 10 3 and ~ when the neutrino goes outwards from the center 
to the sun surface. 

In the general case Eq.(|2|) must be solved numerically ([^, ||). Comparation with 
numerical studies and qualitative analyses show that for slow (in the neutrino time scale 
) and monotonically varying matter density, the adiabatic approximation can be applied 
except, possibly, in some small "resonance" region (MSW effect). This resonance effect is 
essentially the same as the non-adiabatic crossing of energy levels in diatomic molecules 
studied already 60 years ago. The overall transition probability depends essentially of the 
existence of this resonance region and not on the detailed matter density function. The 
aproximate transition probability was computed for two species in two classical papers by 
Landau and Zener || supposing that the density varies linearly in the (small) resonance 
region, its formula, adapted to the neutrino context, is: 



Pee = exp 



Am 2 s 2 29 

-7T 



2EX c26 

The point t res where the resonance appears is given by 

Am 2 cos 26 



(3) 



p(t r 



E 



6 is the 2-neutrino mixing angle. 

Apart from this approximate results, exact solutions have appeared for particular forms 
of the function p: for linear densities, in terms of Weber-Hermite functions ([Q, ^]); for 
functions of the form 

p(t) = C(l + tanh(Ai)) 
in ||, and for exponentially decaying densities 



p(t) = ce 



-\t 



m 



10]. All these functions for the density are of interest in the solar case, the last two 
reproduce rather well the real solar density except in the inner core and to a less degree 
in the surface boundary. All of these solutions are valid for 2 species and all of them are 
obtained reducing Eq.@ to a second order differential equation identifiable with one of 
the classical equations of the mathematical physics. 

In this paper , using a completely different method, summing a standard time-dependent 
perturbative expansion, we will give an exact, analytical solution for an arbitrary num- 
ber of species and for arbitrary (possibly non-unitary in case of sterile neutrinos) mixing 
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matrix using the exponentially decaying function. Using the limit t — ► oo, which is appro- 
priated for the case of detection at Earth, we will give specially simple solutions (this is 



something that papers M, 10] fail to give in the 2-dimensional case). 



2 An exact solution for 5(t) « exp(— t). 

We will get an exact solution for the neutrino evolution equations in the special case where 
the electron density decays exponentially. Standard non-stationary perturbation theory 
will be used; we will give general expressions for the n-th term and sum the full series 
(regardless of mathematical convergence problems, we will see a posteriori that the series 
effectively converges, and the result is physically meaningful). 
The evolution operator of the differential system 



is such that 



obeying the integral equation 



id t v = H(t)u (4) 
u(jt) = U(t,to)v(to) (5) 



l7(Mo) = l-t / H(r)U(T,t )dT (6) 
Jto 

We are interested in the case where the hamiltonian, of finite dimension k , can be 
decomposed in the following way 

H = H° + V(t) = H°+p(t)V (7) 

H° ,the free hamiltonian, hermitic and independent of time, has eigenvalues {E\, . . . , E}.} 
and eigenvectors {| a >,| b > ...}. We can suppose that at least two of the eigenvalues 
are distinct, the completely degenerate case can be solved in a trivial way. p(t) a scalar 
function and V a hermitic matrix with all eigenvalues but one equal to zero. 
Formally, it is possible to solve Eq.(||) by succesive iterations: 

oo 

U(t,t ) = U^(t,t ) + Y,U^(t,t ) (8) 

n=l 

where the order can be taken as: 

C/ (0) (Mo) =exp(-iH°(t-toj) (9) 
and is the well-known integral 

V (n) = ( _ . } n J ^ dTiU 0^ Tn)V{Tn) . . . £/0 (T2) Tl )F( n )tf°(Ti, * ) (10) 
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The domain of integration is defined by 

r = t > r n > . . . > n > t . 
A V with the special form of interest to us can be always reduced to the form 

V = u- 1 Diag(l,0, ... , 0) u (11) 
for some unitary matrix u and its elements are given by 

V V = Yl u ik 81,1811% = u Xj (12) 
k,l 

The following property will be used later valid for any V of this form: 
VijVji = VuVjj (no summation involved) 

First we will compute U in the t — > 00 limit. And the end we will see that in fact 
we can compute it for any finite t using the properties of the evolution operator. For 
simplicity, we will work in the basis of eigenstates of H°. For practical applications (weak 
eigenstates neutrino transition probabilities) we may want to apply a straightforward 
coordinate change U' = v~ 1 Uv later. 

Trought elementary manipulations we get, from Eq.(|l0|), 

< b I U {n) I a >= {-i) n exp(-i(E b t - E a t ))x 

x J d n T exp(iT n w bkl + . . • + iriui fc ( n „ 1 ) a )V r b fci(r n ) . . • V fc(n _i )a (ri) (13) 

kl,...,k(n— 1) ^ 

With Wk\k2 = Eki ~ Ek2- Besides, using the special form of the potential: 

< b I C/ (n) I a >= (-i) n exp(-i{E b t - E a t ))V ba Efci,...,fc(n-i) B kX ... B^^A^^^^ 
Aki,...,k(n-i) = J r d n Texp(iT n w bk i + ... + iTiw k{n _ 1)a )p{T n ) . . . p(n) 

Bk = Vkk 

(14) 

The multiple sums run over all the eigenvalues of the hamiltonian independently. 

We consider an idealized exponentially decaying electron density profile for the sun: 

p(t) = p exp(-At) 

Where po can be positive or negative. Redefining E^, po we can always suppose A = 1. 
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For this kind of density we can use the following equality, 

rco rx n -i _ 

I n (wi, ... ,w n ) = / .../ dx\ . . . dx n exp V w n x n 
Jto Jto - 

(-l) n exp(t En w n) 



W n (w n + W n -i) ...(w n + Wn-! ...+Wi) 

(valid if Sw^O.Vn) (15) 

This equality can be proved by induction, noting that, with the help of the variable 
change 

£1=2/1 + 3/2; X 2 = V2 

the following recurrence relation holds 

I n (wi, . . .,w n ) = h(wi) x I n -i(wi + w 2 ,w 3 , ...,w n ) 

that together with 

-exp(u> t ) 

h(w) = 

w 

proves our result. 
In our case 



w n + w n -i + . . . + w n - j+1 = iw bkl - 1 + iw klk2 - 1 + • . . + iwk(j-i)kj ~ 1 

= iwbkj ~ 3 

so the (n-l)-tensor Aj.j can be factorized 

to) 



(16) 



^l, fc 2,... )fe (n-l) = eiWbat ° S b e a -U AklAk2 ■ ■ ■ (17) 



with 



A k(m) - " ( 18 ) 

*W&fc(m) - ™ 

the fact that this factorization is possible is the key for the solution of the problem. Note 
also the importance of the ordering in the definition of the A's. 
So, 

n-l 

^2 B kl ■ ■ ■ B k (n-l)A k l ■ ■ ■ ^fe(n-l) II X! ^fc(m)-4fc(m) = Il^ m ( 19 ) 

m=l all eigenvalues m 

What remains is the computation of the expression FJ m fm- Such computation is easy 
even in a general case but rather messy. We will do it first for the 2-dimensional case for 
the sake of clarity. 



6 



2.0.1 2-dimensional case 

For a two-dimensional case and taking b=l 



_ B x B 2 _ -1 (-iB lWl2 + m) 

Jin — T - — / • i \ V ZL V 

— m iw\2 — m m {—iwi2 + m) 



and the product of them 



17/ -iBTlMw (21) 

i = V m -B l( n-l)! G9) (B) (21) 



where we have used the Pochammer symbol Eq. (p9|) and defined /3 = —iww Inserting 
this product in the previous expressions and recalling the series expansion of the Confluent 
Hypergeometric function Eq. (68) we can write immediately the value of the diagonal 
elements of U 



<b\U\b>= 

= e-^-o) (l + ^EH)" ( "'° eXP ^ 0)r ff fm) 

\ n=l ' m=l / 

iEb (t-t ) ( 1 + Ytey^ (-ipoexp(-t )) n (B b (3) {n) \ 
V B *tl Ul fa J 

= e -iE b (t-t ) lFl ( Bb /3,/3,z) (22) 

here f3 = —iw ba and z = — ipo exp(— 1 )- 

For the general case a ^ b is neccesary to work a little bit more to give a closed 
expression. In this case it appears an additional factor in each term of the series which 
doesn't allow for the immediate identification with any known function: 



< b | U I a >- 

\n n—l 

I m 



e -i(E b t- Ea t ) ( + e i(E b - Ea)t0Vba ^ ( _- r (-Poexp(-to)) w g f . 



, (3 — n 

n=l ' m=l 



-i(E b t-Eato) c i(E b -Ea)t V ba z " 1 ( B bP)(n) 

B b ^ {n-l)\l3 + n (5 {n) 

Vba 



e -iE h{ t-t )*ba g{Bbj3 ^. z) (23) 

Dh 



as before here (3 = —iw b a- 
The function 



5W + „ G9) (B) (n-1)! (24) 
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is well defined, the series is absolutely convergent for all z; manipulating termwise, and 
using the Eq.([73|) we get the following differential equation for it 



zg' + (3g = az x Fi(l + a(3, 1 + P;z); g(0) = 



(25) 



whose solution is 



g{z) = az' 13 / y\F x (1 + a/3, 1 + 0; y)dy 



(26) 



Using the Expression (71) for the special case c = 7 + 1, we finally write (note that in 
our case the conditions of validity of the previous integral are fullfilled) : 



9i») 



-z iFi(l + a/3,2 + /3;z) 



l + P 

iF 1 (aP,P;z)- 1 F l (aP,l + P;z) 



(27) 



The remarkable second equality can be deduced from the expressions (73-77). So 
< b I U I a >= 

Z 1^1(1+^,2 + ^2) 



-iE b (t-t ) Vba B b 



B b l + P 



< a\U \ b>~- 



e -iEa(t-t )Y^Jk^ z lFl (l-B a p,2-p;z) 



B a l-P 



e -iK(*-to)_^_ z e* lFl (i _ (1 _ £ a )/3, 2 - /3; -*) 
e -iK(t-t )J^_ z B h p,1- P;-z) 



1-/3 



e -iE a (t-t ) Jab_ z e z lFl *(i + Bbj3j2 + P;-z*) 



l-P 



e -iEa(t-t ) _Y?b_ z e * lFl *(i + Bb p j2 + p;z) 



l-P 



(28) 
(29) 



Where we have used the Eq . (|73|) . The last Formula ( P9|) is only valid for z purely 
imaginary, in the case of general complex z (for example given po complex ) Formula (|28| ) 
is valid. 

We have completed so the computation of U. In summary, the elements of the operator 
U in a basis of eigenvectors of H° are 

U(t,t ) = exp-iH (t - t )U red (p ,t ) 



U r ed{PO,to) 



V21 
Vi 



V12 n 



^ exp(z) G* exp z F* 



(30) 
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with the shorthands 

G = g{z) 
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1 + 



1 F 1 (1 + B 1 (3,2 + p;z), F = 1 F 1 (B 1 p, (3; z) 



As it should be, U{t) obeys the free equation: 

id t U(t) = H° U(t) 

For a reverse sign in po (antineutrino case), the matrix U re d is basically the same except 
for the substitution B\ — ► 1 — By. 



Vl2 Q* 

Vn 



U re d(-pO,to) = exp(-z) [y 

y exp(z) G exp z F 



(31) 



where now F stands for ii ? i((l — B\)(3, (3, z), and G changes similarly. 

For H hermitic U is unitary and we get as a byproduct the following non-trivial identity 
for the absolute value of Hyper geometric functions: 



||iFi(a/3,/3;2)|| a + a(l-a) 
or, using the second Equality in 



1 + 13 



LFi(l + aft2 + /3;2; 



1 (a real) (32) 



| x F 1 {a0,P;z) | 2 + 0—^- | i-Fi (a/3,1 + (3\z) - iFi(a/3, /?; z) | 2 = 1 (a real) (33) 
For the particular case of u being the orthogonal 2-dimensional matrix 

ce -so 

^se ce 

the matrix V is 

( C 2 9 -S9C9\ 



-sece s 2 e ) 

and H° the standard neutrino oscillation hamiltonian given by Eq.(|]). We arrive to (to = 
and restoring A), 

U[f = 1 F 1 (^C\^;^) (34) 

rred _ ™ ~W A r, f , iw 21 n 2 Q n , *«>21 , ~ipO 



UlT = TO ^1 + ^^,2+7;^ (35) 

1 + IW21/A \ A A A J 

We are interested in the expression of the evolution operator in the basis of eigenvalues 
of V, it is straigthforward to compute V = uUu^ 1 . After some algebra, applying the 
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identities we have seen before and averaging out time dependent terms, the probability of 
electron permanence is 



F, 



ill 



1 - S 2 9 1 + C29 



/ iw 2 i 



c 2 e, 1 + 



tw 2 i -ipo 



X 



A 



(36) 



An expression equivalent to the previous one has been obtained already by [11] in terms 
of Whittaker functions. However our expression is simpler and more compact because we 
use the infinite time limit. 

In the limit po — ► or free case, z — > and the hypergeometric functions go to 1, U 
becomes as expected: 

U = exp -iH°t 

In the limit of vanishing mixing angle, C 2 6 — > 1, tan# — ► and \F\(z) — > expz in (34), 

so 



1 



In the limit w 2 i/\ = Am 2 /2EX -> 0, xF\ -> 1 in (M) and 



1 



1 



S 2 26 



(37) 



(38) 



In the the limit w 2 i/X >> 1 (adiabatic regime) we distinguish two cases: if w 2 \/X z \ 
we have \F± — > exp(—iC 2 9po/X) and Formula (|38|) also applies; if W21/X ~| 2 | a resonance 
occurs and the probability drops abruptly. Some asymptotic formulas exist for this case 



valid except in the central region of the resonance (see [12 



In the limit of small (but not vanishing) mixing angle, we are left in principle with the 
expression for the electron permanence probability (computed ignoring diagonal terms in 
U) 



;i - 2C 2 es 2 e) 



11 I A ' A ' 



-ipo 



(39) 



In the figures we show a comparison between formulas (£36|) and (|39|). We see that even 
at very small angles Eq.p^) is not a very good approximation of ( |36[ ) as we could expect. 

In Figure ((l|) we plot | \F\ \ 2 (Eq.(^9|)) as a function of its parameters keeping z 
constant. Increasing values of the x-axis implies a larger mass difference or a smaller 
neutrino energy. We observe that the probability is near 1 ( adiabatic evolution ) except 
for a very concrete region of the parameter space (MSW effect or existence of a resonance 
layer). The starting point and width of this region (but not its end point) depends on the 
parameters involved, in particular cos 2 (9. 

Figure @) is equivalent to Fig.((l|), this time we plot the exact expression for the 
e-e probability, Eq. (|36|) , The global behaviour is exactly the same as before, the local 
oscillations have disappeared now however. Figure (^) show the behaviour of the resonance 
region also for large mixing angles. The continuos curve is our exact solution, the dashed 
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curve the Parker approximate formula as given in (ITj]). Figure (||) show the behaviour 
for antineutrinos (reverse sign for pq in z), in this case the resonance is absent. 

In a second set of plots (figs.(|5],||)) we keep constant the parameters and vary z or 
equivalently the production point of the neutrino: larger | z | or po means larger density 
at the creation point which implies that this one is nearer to the center. We see again there 
is a transition between two well defined zones: for a very far creation point, the neutrino 
doesn't pass through any resonance layer, adiabatic regime conditions always apply and 

Pee ~ I- 

In order to study in detail the existence and properties of the resonance region it would 
be convenient to have general formulas for the zeros of the Hypergeometric functions, as 
function not only of its argument but also of its parameters. Unfortunely very little is 
known in a general case. There are some results (see appendix and reference therein) 
about the real zeros of \Fi(a, b; x) with a, b real. 



We note that the expressions (pOI) and ( 39 ) are still formally valid for a purely imaginary 



A (we take for granted that we can take convenient limits 3ft w n — > in the integral (|l~5|)). 
Expression (|3^) however is not valid any more because we used the hermiticity of H (and 
unitarity of U) in computing it. The extrema of | \F\ | 2 are given essentially by the zeros 
of i-Fi(l + a, 1 + b;xj); for A purely imaginary, argument and parameters in iiq become 
real and we can apply the bounds in Appendix A given by the Eqs.(78-3(]) . This bounds 



can be used to deduce the regions in the parameter space which allow for the existence of 
such extrema. 

2.0.2 3- and k-dimensional cases 

For three neutrino species, we can write similarly to Eq.(^0|) 

f _ B i , B 2 , B 3 _ -1 (-oi + m) {-a 2 + m) , . 

— m iwyi — m iw\3 — m m {—iwi2 + m) {—iw\3 + m) 
where oi, 02 are the roots of a certain 2-degree polynomial, they obey: 

ai+a 2 = i (wn(l ~ B 2 ) + wi 3 (l - B 3 )) 

(41) 

aia 2 = -B1W12W13 
and the product of them 

TTf - 1 (- 1 )"' 1 (-«l)(n)(-«2)(n) 

i = y m B X (n-i)\ (0o (B) o%) (B) 1 ; 

with P12 = — 2^12,13. In the simplification of the previous equation has been important 
the expression for a\a,2 in Eq.(^). Note that this last expression is only valid if wyi and 
w\ 3 are both different from zero. 

For the diagonal elements of U, we can write immediately, in terms of the generalized 
Hypergeometric functions (Definition (f70|)) 
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<b\U\b>= e - i£ M*-*o) 2 f : 



2^2(-ai, -a 2 , -iw bj ki, -iw h , k2 ; z) 



(43) 



with 



at + a 2 =i (wb,klO- ~ Bki) + w b ^ k2 (l - B k2 )) 



(44) 



a\a 2 



-B b w bM Wh k2 



k±, k 2 are labels for the two eigenvectors different from b. 

For the < b \ U \ a > (a ^ b), we can repeat exactly the same path as we did in the 
2-dimensional section. We arrive to (we do it already for the general case): 

zg' + Pig = I 1 '" ° fc ~ 1 z 2 F 2 (-ai + 1, .. .,-o fc _i + + 1, . . .,(3 k -i + l',z), g(0) = 



Pi ■ ■ ■ Pk-1 

with Pi = —iw ba ,p k = —iw bk , where k ^ a,b. The solution of which is 



(45) 



g(z) 



«i 



Pi 

a i 



Pi 
«i 



Pi 



^^-z 131 [ dww Pl fc_iF fe _i(-ai + l,...,-a fe _i + l; Pi,..., P k -i,w) 



Pk- 



Pk-i Jo 

Ofc-l z 



z / dss 11 



k-l 



F fc _i(-ai + 1, . . . , -a k -i + l;Pi,...,P k -i;sz) 



Pk-i 1 + Pi 



«-i-Ffe-i(-oi + l, • • • , -Ofc-i + 1; /?i + 2, p 2 + l, . . . , p k -i + 1; z) 

(46) 



Where we have used the expression Eq.(72). for \i = 1 
So 



< 6 I Z7 I a >= e - i{ - Eht - Eat ^V hl 



' 1 + Wfea 



2 F 2 (-ai + 1, -a 2 + 1; -iw b , k i + 2, -iw b)k2 + 1; z) 

(47) 



ai, a 2 fullfill the equations (44) with fci = a, fc 2 are labels for the eigenvectors different to 
b . 

Let's take for u a general unitary 3-dimensional parametrized as follows 



1 

h-io Ci) Sip 

.0 -Si/> Cty. 
The potential V is 



1 

e iS 
e 








C(f) Sep' 
1 
-S(p C<p>, 



Cuj Suj (T 
-5w Cu; 
1. 



The last equality is valid only in our particurlar case (or for any real u). 



(48) 
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u\ a {a = 1 — 3) = (C(pCuj, C(j)Suj, S(p) 
Note that the phase S and the angle ip don't appear in V. 

In the basis where U is diagonal U' = ullu^ 1 . In particular, for the e-e transition we 
have 

An = t/n = E Mu'^kiUik = ]T V ki U k i = trUV 

kl kl 

In the small mixing limit the terms Vki are nearly zero for k, I 7^ 1. Averaging out time 
dependent terms, we write approximately: 

Pee =| U' u | 2 ~ C 2 4>C 2 OJ I U U (= 2F2) | 2 (49) 

with obvious arguments for the hyper geometric function. The computation of a exact, fully 
simplified expression for C7{ 1 requires the algebraic manipulation of several 2 F 2 functions 
and is not of major interest. For the expected degree of approximation of this formula, 
see the comments correspondent to the 2-dimensional case. 

Now we are going to study the behaviour of Uu for two different limits. Let's suppose 
that the first and second eigenvalues are nearly degenerate: wu — ► 0, in this case 

= iw 13 (l-B 3 ) 
= -B1W12W13 —> 

and 

Uu oc l-—( 1 Fi(-ai,-iw 13 ;z)-l) 

= 1 - £1 1013(1 - iF 1 (-iwi 3 (l - B 3 ),-iwi 3 ;z) (51) 

The 3-oscillations has been reduced essentially to a 2 dimensional problem. We can ap- 
ply to this last formula the arguments of previous section for finding different limiting 
behaviours as function of the size of wi 3 . 

It is interesting also to study what happens when one of the eigenvalues is much bigger 
than the others (but none goes to zero) , so some of the to's goes to infinity. Let's suppose 
w\ 3 — > 00. In such a limit , solving previously the second degree algebraic equation, we 
have 

C _ iw 12 (l-B 2 ) 

< 2 (52) 

[ a 2 = iw 13 (l - B 3 ) = iwi 3 (Bi + B 2 ) 



Thus 



lim Uu oc hm 2 F 2 ( ^ 12 ^ 1 — — , iw\ 3 {Bi + B 2 ), iw 12 , iwi 3 ; z) 

= ^ C^-^ ^z^ + B,)) (53) 
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so, as expected, the system behaves effectively as a 2-dimensional system. Note that the 
2-dimensional submatrix of u is not unitary, so the sum B\ + B2 is not equal to unity 
anymore in Eq.(|53|). This procedure is specially useful for the study of oscillations in 
models with extremely heavy sterile neutrinos as those inspired in SO(10) GUT theories. 
For k dimensions, we can write similarly to Eq.(|20|) 



B X B 2 B k (-l) fc ((-g 1 )+m)...((-a fc _ 1 ) + m) 

Jm = + " + • • • + " = r - ; r -j—. ; r- (54 

— m iw\2 — m iwik — m m \—iWx% +m) ... (—iwik + m) 

where ax, ... , afc-i are the roots of the (k-l)-degree polynomial, 

k 

p(m) = ^B s Y[(m-iwxj) (55) 

The coefficient of the power of greatest degree, Co, and the constant term, Ck-i of this 
polynomial are: 

co = E^ = 1 

s 

k 

cfc-i = p(o) =y. b * n(-»>ij = (-i) k ~ iB i n w v 



the product of the roots is then 



ax----- a k _x = (-l)*" 1 ^ = i k - l Bx JJ 

Based in the previous formula and assuming all the Wij different from zero, we can 
generalize the 2 and 3-dimensional expressions for the product, 

TT f 1 (-1)" (-ai)(n)---(-afc-l)(n) , , 

14 Jm Bx(n-1)\ (A)(»)-03*-i)(») 

with j3 s = —iwbs- 

As it can be deduced trivially, the matrix U for k dimensions is the product 

U(t) = exp-iH (t-to)U red (H ,V,po,t ) 

The elements of U re d being essentially Generalized Confluent Hyper geometric Functions 
of order k-1. The diagonal elements are given in a basis of eigenvalues of H° by: 

< a I U red I a >= k-l F k-i(-ax, -a fc _i,/3i, . . .,p k _i;z) (57) 

with the a's and /3's defined as before, and z = —ie~ to ^ x po/X. 
For non-diagonal elements: 
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< b | U red | a >= — - fc _iF fe _i(-ai + l,...,-a fc _i + l;/3i+2,...,/3 fc _i + l;z) (58) 

Vbb 1 + Pi 

Note that for each individual element /3's and a's are differently defined. 

As it was pointed out before the previous expressions are valid only for the completely 
non-degenerate case: all the Wij different from zero. The general case will be treated 
elsewhere, it is easy to see that the order of the hypergemetric functions appearing there 
will decrease in one unit for each of the w's equal to zero. 

The unitarity of U for any hermitian H induces a tower of identities between the 
absolute values of the Hypergeometric Functions generalizing the expressions ( j3^]33| ) found 
for the 2-dimensional case. 



2.0.3 Evolution for finite time. 

Until now we have considered the evolution of the system from a finite to until infinite time. 
Now, using general properties of the evolution operator, we will compute the evolution for 
arbitrary finite intervals. 

Any U, for any intermediate time t\, obeys the composition property: 

U(t,t ) = l/(Mi)tf(fi,*o) 

Taking the limit t — > oo in the above expression we arrive (taking to = for simplicity) 

U oo (t,0) = U oo (t,h)U(t 1 ,0) 

By the subindex oo we want to remark we have computed U as in the previous sections, 
using some limiting process. So , as U has an inverse, 

u(t,o) = [/-V.^ocGM) 

= uljto = ty Ho ^e-^u red {tv = o) 

= Uljto = t)e~ iHot U red (t = 0) (59) 

Where we have effectuated some relabing to let clear the time dependence. \i , appearing 
at first as a new arbitrary parameter, disappears from the last expression. U red is the 
matrix depending on the diverse parameters of the system, including to (compare with 
Eq.©). 

The matrix defined by Eq. ( |59| ) satisfies by construction the initial condition U(0, 0) = 
1. In addition we must check it satisfies the Equation (Q): inserting the last expression we 
obtained, and eliminating the common factor U red (to = 0), we get 
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- id t e lHot U red 



= e iHot U red (H + p(t)V) 



e iHot H U red - ie iHot d t U red = e iHot U red (H + p{t)V) 

[H ,U red ] = id t U red + p(t)U red V (60) 

the dependence of U red on t comes from the substitution to = t. 
If 

[H ,U red (t)] »0 (61) 

for example if the difference of the eigenvalues of Hq is relatively small, then Eq. ( |60[ ) has 
the trivial solution: 

Ured{t) = exp-iV [ g(p)dp (62) 
Jo 

and the U matrix becomes 

U(t, 0) = exp —iH(jt U red {t) = exp —iH$t exp — iV \ g{p)d\x (63) 

J o 

this is a highly non-trivial expression valid for arbitrary Hq, V and git). Some pertubation 
approximation expansion in AE can be implemented for better approximations. 

In our particular case it is easy to check that the previous Eq.(|60[) is really fullnHed. 
It converts to a system of algebraic- functional relations between Hypergeometric func- 
tions. For simplicity we take the 2-dimensional case. Eq. (^0|) is equivalent to only two 
independent relations, with the same notation as in previous sections: 



d z F = F V u + | V 12 | 2 G 
F = (- Vl2 + (3/z)G + d z G (64) 

with the help of Eq.(]73|), we obtain instead two algebraic equations that can be proved 
to be identities with Eqs.(|77|). We just prove in detail the first one. First we identify 
Vxi = a, = a(l — a), so the first term is 

d z F(af3,f3,z) = a]Fi(a/3 + 1,(3 + 1; z) (65) 

the second term is 

F Vn+ I V12 | 2 G = aiFi(aP,P;z) + (l-a)(iF 1 (af3,P;z)-xF 1 (ap,p + l;z)) 

= 1 F 1 (aP,/3;z)-(l-a) 1 F 1 (aP,P + l;z) (66) 
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In the first line we have used Eq. (|33[). The last expression is just the identity fl7q) . 
In summary, we have proved that the evolution operator for any finite time is 

U(t,0) = ul d (t)e- iHot U red (0) (67) 

with U re d given in the previous sections. Expressions can be given now for the transition 
probabilities in the mass or weak basis as we did before. However the expressions are 
rather involved and not so illuminating. For almost any practical purpose we can stick to 
the expressions given in the t — > oo limit. 

3 Conclusions and further discussion. 

Exact expressions for the transition probabilities of neutrinos propagating radially in the 
sun are given. These are valid for an arbitrary number of neutrino species. The solution, 
very compact in terms of Hypergeometric Functions, is very suitable for analytical studies 
of resonance and for systematic numerical approximations. The computational task of 
computing numerically this kind of functions makes however very impractical the actual 
use of them in concrete applications for the moment. 

The perturbative expansion method used, that in this particular case gives an exact 
solution, is well suited to produce approximate results for the same matrix form of the 
potential (all eigenvalues but one equal to zero) but for arbitrary form of p(t). The 
consideration of a finite number of terms in the expansion is expected not to be useful 
but yes the summation for all orders of appropriate leading terms. Another possibility is 
to consider the difference between an arbitrary p(t) and our exponential as small and do 
perturbation theory around our exact solution. In principle we expect this procedure to 
be convenient at far distances from the creation point, as we have seen the probability 
transition depends more on the existence of a resonant region than on the detailed shape of 
the potential. This would make possible to compute formulas for non-radial propagation 
and for the regions where the exponentially decaying density lose validity: central core 
and solar surface. 

The phenomenon of resonance (or MSW effect) is remarkable in itself. The present 
solution offers a starting point for its analytical description. We note the similarity of this 
effect with the original form of the Anderson Localization effect presented in |15|] . This 
similarity is specially evident under the treatment done in this work. In both cases the 
locking of a system in one of its quantum states except for a certain range of the param- 
eter space arise as a property of the solution of a coupled system of ordinary differential 
equations. As in its case (which physically correspond to the evolution of a particular 
site in a random lattice) the proper energies (or its differences) of each mode play an 
important role, however here these are not needed to be stochatically distributed. This 
condition of randomness is essential in |15[| , the results in this work induce us to think 
whether such randomness is a mere mathematical convenience more than having a deep 
physical meaning. 
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Apart from physical interest, our solution has intrinsic mathematical interest in the 
Theory of Differential equations and of Special Functions. A general theory of equations 
of the form 

d t x = {A + f(t)B)x 

with A,B constant matrix is missing in spite of the fact these systems are the "next-step" 
in complexity from all-constant coefficient equations. Only solutions for a few examples 
of these equations are known. 

The generalized Hypergeometric functions are shown to be the asymtotic solutions of 
a simple difererential equation in a systematic way. New identities for the absolute values 
of these functions are derived. The generality of the results obtained here induces us to 
consider the utility of defining Hypergeometric functions of matrix parameters \Fi{A, B; z) 
with z complex, similar extensions exist already: generalizations with complex parameters 
but matrix argument or parameters and argument defined in arbitrary finite fields for 
example. Let's note finally that making the change 

y = exp(-i) 

our system becomes of the form 

d y x = ^— + B^j x 

so we have found also a solution for this system, the Hypergeometric functions appear 
now as regular solutions for y — ► 0. 

A Appendix: some formulas about Hypergeometric Func- 
tions 



Sec [121 and references therein for all of these definitions and formulas. The Confluent 



Hypergeometric function is defined by 

-S&s (68) 

where the Pochammer symbol is 

{z) {n) =T(n + z)/T(z) (69) 
The generalized Hypergeometric function is defined by 

77i / r t \ f (Ql)(n)---(Op)(n) Z n 

p F q (ai, ...,a p ,bi,...,b q ;z) = 77-7 77^ r (70) 

The integral 
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f x^\t-x) c -^\F l {a,r,x)dx = t c - l V{l) lf. 7) iFi(a,c;t); [»c, 7 >0] (71) 
JO 1(c) 

is a special case of 



o 



(1 — x)^ 1 - ) x bl 1 pF g (a 1 , . . . , a p ; b\, . . . , b p ; ax)dx 



i-pF g (ai, • • • , a p ; /i + 6i, 6 2 , • • • , a); 



r( M + 6i)' 

(»0i,6i)>0, p<g+l) (72) 
Some other important formulas are: 

d i-Fi(a, 7; z) a „ . , . , 

\ = ~ 1F1 l + a,l+ 7 ;z 73 
az 7 

1 F 1 (a,~f;z) = expz 1^1(7 - a,7; -z) (74) 

-lFi(a + l,7 + l;z) = iF^l + a, 7; z) - iFi(a, 7; z) (75) 
7 

aiFi(a + 1,7 + 1; z) = (a - 7)1^1 (a, 7 + 1; z) + 71^1(0,7; z) (76) 

aiFi (a + 1,7; z) = (z + 2a - 7)iFi(a,7;z) + (7 - a)iFi(a - 1, 7; z) (77) 



From [16| we know that the real zeros Xj of iFi(a, c; x) for a,c real satisfy the bounds 



(c - 2a) - 2y/(a(a - c) - c) < Zj < (c - 2a) + 2^/(a(a - c) - c) (78) 
The smallest real zero x\ = x m in satisfy 

c(c + 2) 

Xmin ^ ^ I'** J 

c — la 

Appliying the same bound to the expression ([74]) we deduce the lower bound for the 
maximal real zero x max 

[x min <) — < (80) 

c — 2a 
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Figure 1: The function | iFi(Bi(3,i(3, — i550.2) | 2 for two different B's near 1, as a function 
of (3. This corresponds approximately to the oscillation probability for a neutrino produced 
at r/r = 0.25 as function of j3 = Am 2 /2EX. 
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Figure 2: P ee (two neutrino species, Formula 36) for a neutrino produced at r/ro = 0.25 as 
a function of f3 = Am? /2EX and two different mixing angles: B(= cos 2 9) = 0.99, 0.999. 
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{B=0.50} 



{B=0 .70} 




Figure 3: As in Figure (g) but for a neutrino created near the surface (| z |« 20) as 
a function of the mixing angle. The dashed line corresponds to the Parker approximate 
solution (as it is given by Formula 2.18 in Q). 
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(B=0.70) 




10. 100. 



10. 100. 



Figure 4: As Fig. (||) but for antineutrinos (reverse sign in the argument z of the hyper- 
geometric function). In this case there is not resonance (it is possible to distinguish some 
"anti"-resonace behaviour for larger mxing angles). 
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Figure 5: P ee (Eq. (|36|)), as a function of the creation point (in solar radius fraction) for 
a fixed mixing angle cos 2 9 = 0.99 and two different (3 = Am 2 /2EX = 1,20 (respectively 
upper and lower curves ). 
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Figure 6: Function | iFi(iB/3,i(3,z) | 2 as a function of | z |, for different values of the 
parameters. B=0.99, the upper and lower curves correspond respectively to (3 = 1 and 20. 
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